A solvable non-conservative model of Self-Organized Criticality 
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We present the first solvable non-conservative sandpile-like critical model of Self-Organized Critical- 
ity (SOC), and thereby substantiate the suggestion by Vespignani and Zapperi [A. Vespignani and 
S. Zapperi, Phys. Rev. E 57, 6345 (f998)] that a lack of conservation in the microscopic dynamics 
of an SOC-model can be compensated by introducing an external drive and thereby re-establishing 
criticality. The model shown is critical for all values of the conservation parameter. The analytical 
derivation follows the lines of Broker and Grassberger [H.-M. Broker and P. Grassberger, Phys. Rev. 
E 56, 3944 (1997)] and is supported by numerical simulation. In the limit of vanishing conservation 
the Random Neighbor Forest Fire Model (R-FFM) is recovered. 



PACS numbers: 64.60.Ht, 05.65. +b, 02.50.-r 

The role of conservation in SOC-models is an old is- 
sue and is still unsettled. The number of non- 
conservative models which are definitely critical is, how- 
ever, strikingly small. The Random Neighbor Forest Fire 
Model (R-FFM) is one of them, while the Random 
Neighbor Olami-Feder-Christensen model (R-OFC) has 
been shown not to be critical in the non-conservative 
regime |>KJ. The nearest neighbor OFC model is com- 
monly accepted to be critical in the conservative limit, 
but whether this model is critical in the non-conservative 
regime is still debated 

In J8|,|| it has been suggested that non-conservation in 
the microscopic dynamics can be compensated by an ex- 
ternal drive in order to re-establish criticality. Applying 
this concept directly to a model known to be non-critical 
in its original definition provides the ideal basis to iden- 
tify the effect of such an external drive. In this letter such 
a model is defined and solved semi-analytically. The re- 
sults are compared to simulations and the (trivial) critical 
exponents are extracted. Several limits are discussed. 

The model, which is derived from the FFM |l(| and 
the Zhang model jy], has three main parameters: N is 
the total number of sites, which diverges in the ther- 
modynamic limit. The number of randomly chosen 
"neighbors" is given by n, where n — 4 in all exam- 
ples, corresponding to a two dimensional square lattice. 
The conservation parameter is a. The degree of non- 
conservation is then 1 — not, as it is shown below. Each 
site i £ {1,2, ••■,7V} has associated a value Zj for its 
"energy" . Sites with < Z{ < 1 — a are said to be "sta- 
ble" , sites with 1 — a < Zi < 1 are called "critical" and 
sites with 1 < Zi are "active" . Negative energies are not 
allowed. The probability density function (PDF) for the 
variable Zi is P(z) with z £ [0, 1 [ and is defined only when 
no sites are active. P(z) contains most of the stationary 
properties of the model. 

The dynamics of the model are defined as follows: Af- 
ter an initial choice of zi (i = 1, 2, TV) from a uniform 
distribution in the interval [0, 1[, the model is updated 
by repeatedly (i) "driving" , (ii) "triggering" and (iii) "re- 



laxing" the system. During the drive (i), i — 1, (1/0) 
sites are chosen randomly ((1/0) = p/f in the notation 
of 1 10 1 ) , one after the other, and their energy Zi are set to 
1— a, if the site is stable, otherwise Zi remains unchanged. 
Subsequently one random site j is chosen and if it is crit- 
ical, the system is triggered (ii) by setting the energy of 
the chosen site to 1, i.e. making it active (initial seed). 
Otherwise Zj remains unchanged and the model is driven 
again by repeating (i). As long as N and (1/0) are finite, 
the system will escape from the driving loop sooner or 
later. In the thermodynamic limit this is ensured by a 
non vanishing density of critical sites. 

During the relaxation (iii) the energy of each active site 
i is redistributed according to the conservation parameter 
a to n randomly chosen sites j and the energy of Zi is 
then set to 0: 



Zj + azi 







. j ,j ., \j (1) 

Each visit or "toppling" (Q) defines a microscopic time 
step and dissipates exactly (1 — na)zi energy units. The 
sites j are chosen randomly one after the other and are 
not necessarily different. In the thermodynamic limit the 
probability of choosing a target site which is already ac- 
tive or was already charged during the same avalanche, 
vanishes and therefore the order of these visits is irrel- 
evant. In this very restricted sense the model might 
be considered as "Abelian" . In contrast, sites in finite 
systems have always a finite probability to get charged 
more than once. Nevertheless, this probability decreases 
rapidly with increasing system size. 

The number of active sites relaxed by (Q) defines the 
avalanche size s, which is always positive due to the initial 
seed. In the stationary state the avalanches must dissi- 
pate, on average, the same amount of energy as is sup- 
plied by the external drive and the initial seed. The av- 
erage dissipation depends on the avalanche size weighted 
average energy of active sites z act , which is equivalent to 
the average energy of active sites per toppling. Therefore 

(1 - na)-z— t ( S ) = (1/0)^(1 - a - ~t) + (1 - X) (2) 

Pc 
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must hold exactly in the stationary state even for finite 
systems and does not introduce any approximation. Here 
(s) is the average avalanche size, p^t (pH) is the density 
of stable (critical) sites (the drive stage is, on average, 
repeated 1/pT times), zit and z^ are the average energy 
of stable and critical sites respectively. As in || the only 
crucial assumption is that (s)/N as well as (1/9) /N van- 
ishes in the thermodynamic limit, which turns out to be 
entirely consistent with the results. This assumption al- 
lows us, for example, to assume the distribution P(z) to 
be essentially unaffected by external drive or relaxation 
for sufficiently large systems. 

From (||) it is clear that in general (s) diverges for 
diverging (1/9) or vanishing dissipation rate 1—na. From 
the microscopic dynamics it is clear that there is always 
a non-vanishing fraction of sites with z — 0, therefore 
(1 — a) — z^ is finite and a divergence of (1/9) entails a 
divergence of (s), which is a sign of criticality. 

In the following outline of the actual calculation, which 
is adapted from the PDF's of the model arc derived. 

After an avalanche, each site belongs to one of m + 2 
classes, where m — [l/a\. The index k — 0,1,..., m of 
the class indicates the number of charges received from 
other toppling sites since their last toppling, while k = 
m + 1 indicates the class of sites, whose energy has been 
set by external drive. A site charged more than m times 
must be active. For each of these classes a conditional 
distribution function Q k (z) is introduced, describing the 
distribution of energy among non- active sites, which have 
been charged k = 0, m times or externally driven, k = 
m + 1. The distribution of sites which have not been 
charged since their last toppling, Qo(z), is a delta peak 
at z = 0. For convenience the normalization of Qo(z) is 
chosen to be unity and all other distribution functions are 
normalized relative to class 0. The distribution of sites, 
which have been driven externally and have not changed 
since then, Q m+ i(z), is obviously a delta peak at 1 — a. 
If the fraction of externally driven sites is g, P(z) can be 
written as 



m+l 

P(z)=NY,0-9)Qk[z) 

k=0 



(3) 



with appropriately chosen normalization M < 1. The up- 
per bound for the energy of an active site is the geomet- 
ric sum 1 + a(l + a(- ■ ■)) — 1/(1 — a), neglecting double 
charges. Therefore, the support of the distribution func- 
tion of active sites C(z) is [1, 1/(1 — a)[. If this distribu- 
tion is normalized, the expected increase per avalanche 
in the class k > is given, in the thermodynamic limit 
(where multiple toppling can be neglected), by the con- 
volution 

-l/(l-a) 

dz'C(z r )Q k - 1 (z-az') , 



n(s) 



(4) 



where the factor n(s) takes into account the expected to- 
tal number of charges. There are three different ways in 



which the classes k < m + 1 may be decreased: 

1) By charges, Q k (z)n(s) 

2) By external drive, Q k (z)9<((l - a) - z)(l/0)^~\ 
where 9 < is the Heaviside unit function with 6 I< (0) = 0. 

3) By initial seed, Q k (z)9> (z-(l-a))/p^, where 0>(O) = 
1 correspondingly. 

Adding these contributions together and assuming sta- 
tionarity leads to m equations for Q kl k = 1, • • • , m: 



Qk(z)l(z) 



where 



V(i-") 

dz'C(z')Q k . 1 (z-az') 



(1/9) «<, 



(5) 



l(z) = l + ±Yr9 <((l-a)-z) + 
p c n(s) 

'(*-(!-«)) 



p c n(s) 



(6) 



has been used. For diverging (1/9), the last term in (|^) 
becomes irrelevant and the RHS of (||) is dominated by 
the first term, meaning that the initial seed becomes ir- 
relevant compared to the external drive. It is reasonable 
to restrict the range of a so that single charges cannot 
activate a site, a/(l — a) < 1 a < 1/2 (due to na < 1, 
this is a restriction only for n = 1). Then one of the Eqs. 
in (||) can be written as 



(7) 



due to the particularly simple form of Qo(z). 

Since Qq(z) — S(z) by definition and Af(l — 
g)Q m +i(z) = g5((l — a) — z) from (^), one further equa- 
tion is necessary in order to find m + 3 distributions Q k , 
k = 0, ■ • • ,m + 1 and C(z). This equation concerns the 
construction of the distribution of active sites C(z). Since 
active sites are created due to charging or as the initial 
seed the average distribution of the number of those sites 
per avalanche is given by 

,l/(l-a) 

(s)C(z) = n(s) J dz'C(z')P(z-az') + 6(z-l) (8) 

where the 5-function represents the initial seed. 

Although it is a priori unknown whether there ex- 
ists a stable solution, or whether it is unique, the set of 
equations given above is enough to start an iteration pro- 
cedure in order to find a solution. The scalar parameters 
required are (s) from (g), p^, z^t, ~z^, which are eas- 
ily derived from (||) and z ac t, the first moment of C(z). 
While n and (1/9) parameterize the problem, g remains 
the only unknown quantity, which is found to be 



<J 



pVt(l/9) 
n{s)% + 1 



by comparing the in- and outflow of class m 
externally driven sites, per avalanche. 



(9) 
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All the equations above can alternatively be derived 
directly from the fundamental dynamics. This ensures 
that the solution is exact in the thermodynamic limit 
given the stationarity assumption. 

The implementation of the iteration procedure is 
straight-forward. As a criterion for termination, one 
could check whether C(z), as defined by (||), is prop- 
erly normalized 0], as it can be proven that it must be 
correctly normalized if it is a solution. However, it would 
be sufficient to assume C(z) proportional to the RHS in 
(||). Moreover, in the numerical procedure the quality of 
the normalization of C(z) depends strongly on the reso- 
lution of the grid chosen, whenever C(z) changes rapidly 
as function of z. Therefore convergence of the iteration 
procedure is better verified by checking whether C(z) ap- 
proaches a fixed point, i.e. is invariant under (||). Since 
the distribution is expected to be highly non-analytic - 
there are at least two <5-f unctions in P(z) - sophisticated 
integration routines are inappropriate. For n — 4 the 
procedure quickly converges for 0.07 < a < 0.24, all non- 
pathological initial values tested lead to the same stable 
solutions. Only for small values of (s), when the <5-peak 
of the initial seed starts to propagate through the distri- 
bution, a large grid is required for sufficient resolution. 
The same problem appears close to the commensurable 
limits mentioned below. 
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FIG. 1. Distribution of energy, P(z), for n = 4 and differ- 
ent values of a. Since there is only a <5-peak in < z < 0.2, re- 
sults for z < 0.2 are cut off. Continuous lines indicate results 
from theory (grid size 32000, integrated in 125 bins), points 
represent results from numerical simulations with TV = 10 6 , 
10 6 avalanches for equilibration and 5 10 6 avalanches for 
statistics (125 bins for the histogram). 



In Fig. [j] numerical simulations of the model are com- 
pared to the numerical solution of the analytical ap- 
proach. Although the PDF is very structured, discrep- 
ancies are small and can be reduced by increasing the 
resolution of the underlying grid. 

The distribution C(z) collapses to a (^-function in at 
least two limits. Firstly, when a — > 0, the FFM-limit, 
the probability for a site to become critical due to a num- 
ber of charges vanishes as g^ 1 "")/", where q < 1 is the 
product of the probability that a site receives a charge 
from a relaxing site and the probability that a site is not 
driven externally between two hits. Hence, for a — ► 



the mechanism of "growing by charges" becomes negligi- 
ble and the external drive becomes the dominating source 
for critical sites. The dynamics now become equivalent 
to the FFM: stable sites = empty sites, critical sites = 
trees, and active sites = fires. Furthermore, as a — > 0, 
the support of C(z) becomes smaller and smaller and the 
distribution of active sites is strongly peaked at z = 1, 
collapsing to a ((-function. Therefore, the distributions 
Qk, k > are less smeared out, as shown in Fig. |^(a) for 
a small value of a. Assuming C(z) to be a (^-function, 
one can easily reconcile the results in Q (Eqs. (3) and 
(7)). The assumption that P(z) is unaffected by single 
avalanches corresponds to p, f — ► in the SOC-limit of 
the FFM §J]. 
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FIG. 2. P(z) as in Fig. 
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In the second limit, a — ► 1/n, the model becomes con- 
servative, but more important, a becomes more and more 
commensurable in the sense that a site charged m = n 
times is almost always active and therefore the support of 
Q m vanishes, as it is squeezed between na and 1. When 
a is very close to 1/n, most of the active sites are pro- 
vided by Qm-i and their average energy is just above 1, 
i.e. C(z) becomes more and more S shaped and so do the 
Qi, as shown in Fig. ||(b). The same behavior is obtained 
whenever ka — 1 for k £ N. 

The critical exponents of the model, can be obtained by 
mapping it on to a branching process [p"3[ in order to iden- 
tify the critical exponent 6 = 2, where V(t) oc t~~ h is the 
exponent of avalanche duration. The exponent r = 3/2, 
found by mapping the model on to a random walker along 
an absorbing barrier, is the exponent of avalanche sizes, 
V{s) oc s~ r . Formally these exponents arise only for di- 
verging cutoffs in the distributions, which are controlled 
by the average number of active sites produced per single 
toppling, the branching ratio a. The cutoffs diverge for 
<r->l. 

However, the mapping is non-trivial, except when C{z) 
is a ^-function. This is because a distribution of active 
sites entails a distribution in the branching ratio, i.e. the 
branching probability itself becomes a random variable. 
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This problem is removed by considering the ensemble av- 
erage with an annealed rather than quenched disorder in 
the branching probability, i.e. writing the probabilities 
for one given node (active site) to branch into k new 
nodes (active sites) as 

v(i^k) = {(fyp k (i- P ) n - k ) p (10) 

where p denotes the branching probability (which is 
a function of the energy of the site) and (} p denotes 
the weighted average over the probabilities. Therefore 
a = J2 k kV(l — » k). This branching process is charac- 
terized by the same generating functions as the standard 
branching process [13|] , which becomes critical for a = 1. 
Hence the condition for criticality is 

E*(uK (1 " p)T, ~*) =^)p = i (u) 

fe=0 ^ ' p 

which is the (average) branching ratio, according to (||) 
given by 

r l/(l-a) f l/(l-a) i 

n{p) p =n dz dz' C(z')P(z - az') = 1 - — . 

Ji Ji vV 

(12) 

Moreover C{z) is time or generation dependent, since 
it evolves from an initial distribution which is only a 5- 
function at z = 1. The deviation of C(z) from the final 
distribution decays exponentially fast, which can be seen 
by investigating the Markov chain of the repetitive convo- 
lution of a now time dependent C{z) with P(z) as in (^). 
Therefore the cutoff, introduced by the deviation of a 
from 1, is dominated by the asymptotic iteratively stable 
limit of C (z) only. Since the asymptote is approached ex- 
ponentially fast the transient cannot influence the value 
of the exponents. The same arguments apply for the ran- 
dom walker approach, therefore 6 = 2 and r = 3/2 is true 
for all a <e]0, l/n[. 

The calculations above are a priori valid only in the 
thermodynamic limit. However, a simulation of the 
model must consider a finite system. Moreover the model 
relies on several assumptions, which entail certain finite 
size scaling: (1/6)n/N (the index indicates the value to 
be measured in a system of size N) as well as (s)n/N 
must vanish for diverging N, while (1/0)n/{s)n must 
remain constant. It is a well known problem in the FFM 
that the number of trees grown between two igni- 
tions is a parameter, (1/9)n, which needs to be tuned 
according to the system size; it is supposed to diverge, 
but its value is restricted by system size. An inappro- 
priatly chosen parameter produces a small value of the 
cutoff or a bump in the distribution function of avalanche 
sizes. Nevertheless, P{z) depends only weakly on (1/6) jy. 
As a more quantitative measure for the "right choice of 



(1/0) jv" w e compared gN to g (see Eq. (0)) in the ther- 
modynamic limit. Assuming a cutoff of order O(N) in 
V(s) of a finite system, the scaling is (s)n = / dsV(s)s € 
OiN 1 / 2 ) and thus (1/0) w G 0(N 1 / 2 ). For a more quanti- 
tative picture one can map the avalanche on to a random 
walker along an absorbing barrier with time dependent 
walking probability p(t) (in the sense of (ll|] a drinking 
rather than a drunken random walker). However, the 
result is comparatively involved and gives only rough es- 
timates for the avalanche size as a function of N and 
(1/0)- 

In summary, a solvable SOC model, critical in the en- 
tire regime of the conservation parameter, has been de- 
fined and the main properties have been derived. The 
critical exponents are as expected the trivial exponents 
of a critical branching process and a random walker. The 
model clarifies the role of the external drive and repre- 
sents an explicit example of the recovery of criticality by 
introducing an external drive. 
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